Correcting for Scattered Light in WFC3/IR Exposures: Manually Subtracting Bad Reads

Learning Goals

This notebook presents one of two available methods to correct for a time variable background (TVB) due to scattered light from observing close to the Earth's limb. This method illustrates how to manually subtract any bad reads from the final exposure read of the WFC3/IR IMA data.

By the end of this tutorial, you will:

• Compute and plot the difference between IMA reads to identify those affected by TVB.
• Correct a single exposure in which the first few reads are affected by scattered light by subtracting those "bad" reads from the final IMA read.
• Compare the original FLT to the reprocessed FLT image.
• Compare the original DRZ to the reprocessed DRZ image.

A second method for correcting TVB (masking bad reads in the RAW image and re-running calwf3) can be found in the notebook Correcting for Scattered Light in WFC3/IR Exposures: Using calwf3 to Mask Bad Reads (O'Connor 2023). This method performs the 'ramp fitting' step in calwf3 and therefore removes cosmic rays from the final image. It only works, however, when TVB occurs at the beginning or end of the IR exposure and may leave sky residuals in regions flagged as IR 'blobs'.

Please note that the FLT products in this notebook are really 'corrected IMA' files and therefore do not include the 'ramp fitting' step in calwf3. The final images will therefore still contain cosmic rays, and these artifacts may be removed using software such as AstroDrizzle when combining multiple exposures.

Table of Contents

Introduction
1. Imports
2. Downloading Data
3. Identifying Reads with Earth-limb Scattered Light
4. Querying CRDS for the Reference File

5. Correcting the Data

• 5.1 Excising "Bad" Reads
  
• 5.2 Comparing FLT Products
  

6. Drizzling Nominal and Reprocessed FLT Products
7. Conclusions
Additional Resources
About this Notebook
Citations

Introduction

WFC3 images can be affected by many different kinds of unwanted, non-science features, such as satellite trails, peculiar cosmic rays, Earth-limb scattered light, reflections off the filter wheel, and more. We call these features "anomalies". The majority of these anomalies are well studied, correctable, and easily identifiable, as discussed in WFC3 ISR 2017-22.

In this notebook, we will walk through the process of correcting WFC3/IR images for anomalies that appear in only some of the reads (see section 7.7 of the WFC3 Instrument Handboook for a discussion of the WFC3/IR MULTIACCUM observing mode). We examine an observation affected by a strong TVB due to Earth limb scattered light affecting the first few reads (section 7.10 of the WFC3 Data Handbook). We will manually subtract the affected reads from the final read of the IMA and copy the updated science and error arrays into the FLT science and error arrays. This method is based on section 5 of WFC3 ISR 2016-16: 'Excising Individual Reads'. This same method can be used to remove individual reads that contain other anomalies, such as satellite trails.

This new FLT will have a reduced total exposure time (and signal-to-noise) given the rejection of some number of reads but will now have a flat background and an improved ramp fit.

Please see the notebook WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023) for a walkthrough of how to identify a TVB due to scattered light.

1. Imports

This notebook assumes you have created the virtual environment in WFC3 Library's installation instructions.

We import:

• os for setting environment variables
• glob for finding lists of files
• shutil for managing directories
• numpy for handling array functions
• matplotlib.pyplot for plotting data
• astropy.io fits for accessing FITS files
• astroquery.mast Observations for downloading data
• wfc3tools calwf3 for calibrating WFC3 data
• ginga for finding min/max outlier pixels
• stwcs for updating World Coordinate System images
• drizzlepac for combining images with AstroDrizzle

We import the following modules:

• ima_visualization_and_differencing to take the difference between reads, plot the ramp, and visualize the difference in images

import os 
import glob
import shutil 
import numpy as np
from matplotlib import pyplot as plt
from astropy.io import fits
from astroquery.mast import Observations
from wfc3tools import calwf3
from ginga.util.zscale import zscale
from stwcs import updatewcs
from drizzlepac import astrodrizzle

import ima_visualization_and_differencing as diff

%matplotlib inline

The following task in the stsci.skypac package can be run with TEAL:
                                    skymatch                                    
The following tasks in the drizzlepac package can be run with TEAL:
    astrodrizzle       config_testbed      imagefindpars           mapreg       
       photeq            pixreplace           pixtopix            pixtosky      
  refimagefindpars       resetbits          runastrodriz          skytopix      
     tweakback            tweakreg           updatenpol

2. Downloading Data

The following commands query MAST for the necessary data products and download them to the current directory. Here we obtain WFC3/IR observations from HST Frontier Fields program 14037, Visit BB. We specifically want the observation "icqtbbbxq", as it is strongly affected by Earth limb scattered light. The data products requested are the RAW, IMA, and FLT files.

Warning: this cell may take a few minutes to complete.

OBS_ID = 'ICQTBB020'
data_list = Observations.query_criteria(obs_id=OBS_ID)

file_id = "icqtbbbxq"
Observations.download_products(data_list['obsid'], project = 'CALWF3', obs_id = file_id, 
                mrp_only = False, productSubGroupDescription = ['RAW','IMA','FLT'])

Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/icqtbbbxq_ima.fits to ./mastDownload/HST/icqtbbbxq/icqtbbbxq_ima.fits ... [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/icqtbbbxq_flt.fits to ./mastDownload/HST/icqtbbbxq/icqtbbbxq_flt.fits ... [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/icqtbbbxq_raw.fits to ./mastDownload/HST/icqtbbbxq/icqtbbbxq_raw.fits ... [Done]

Table length=3

Local Path	Status	Message	URL
str47	str8	object	object
./mastDownload/HST/icqtbbbxq/icqtbbbxq_ima.fits	COMPLETE	None	None
./mastDownload/HST/icqtbbbxq/icqtbbbxq_flt.fits	COMPLETE	None	None
./mastDownload/HST/icqtbbbxq/icqtbbbxq_raw.fits	COMPLETE	None	None

Now, we will copy our RAW file into our working directory to use in this tutorial. We copy the IMA and FLT files to a new directory called "orig/" for later use.

if not os.path.exists('orig/'):
    os.mkdir('orig/')

shutil.copy(f'mastDownload/HST/{file_id}/{file_id}_ima.fits',f'orig/{file_id}_ima.fits')
shutil.copy(f'mastDownload/HST/{file_id}/{file_id}_flt.fits',f'orig/{file_id}_flt.fits')  

raw_file = f'mastDownload/HST/{file_id}/{file_id}_raw.fits'

shutil.copy(raw_file, f'{file_id}_raw.fits')

'icqtbbbxq_raw.fits'

3. Identifying Reads with Earth-limb Scattered Light

In this section, we show how to identify the reads impacted by the scattered light by examining the difference in count rate between reads. This section was taken from the WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023) notebook, which includes a more comprehensive walkthrough of identifying time variable background.

Here we implement a technique to examine the count rate difference between consecutive reads. In this case, we first convert from count rate (electrons/second) back to counts (electrons) before taking the difference, as shown in equation 3 from WFC3 ISR 2018-05.

We compare sky values in different regions of the detector (left side, right side, and full frame). If you would like to specify your own regions for the left and right sides of your image, you can change the lhs_region and rhs_region parameters. Each region must be specified as a dictionary including the four "corners" (x0, x1, y0, and y1) of the region you would like to select. You may want to avoid the edges of the detector which have a large number of bad pixels and higher flat field errors.

fig = plt.figure(figsize = (20, 10))

ima_filepath = f'orig/{file_id}_ima.fits'

path, filename = os.path.split(ima_filepath)
cube, integ_time = diff.read_wfc3(ima_filepath)


lhs_region = {"x0":50,"x1":250,"y0":100,"y1":900}
rhs_region = {"x0":700,"x1":900,"y0":100,"y1":900}

#Please use a limit that makes sense for your own data, when running your images through this notebook.
cube[np.abs(cube) > 3] = np.nan

diff_cube = diff.compute_diff_imas(cube, integ_time, diff_method = "instantaneous")
median_diff_fullframe, median_diff_lhs, median_diff_rhs = diff.get_median_fullframe_lhs_rhs(diff_cube, lhs_region = lhs_region, rhs_region = rhs_region)

plt.rc('xtick', labelsize = 20) 
plt.rc('ytick', labelsize = 20) 
plt.rcParams.update({'font.size': 30})
plt.rcParams.update({'lines.markersize': 15})

diff.plot_ramp(ima_filepath, integ_time, median_diff_fullframe, median_diff_lhs, median_diff_rhs)

plt.ylim(0.5,2.5)
_ = plt.title(filename)

Here, we utilize a few functions from our module ima_visualization_and_differencing. We use read_wfc3 to grab the IMA data from all reads and corresponding integration times. We also implement upper and lower limits on our pixel values to exclude sources when plotting our ramp. We take the instantaneous difference using compute_diff_imas, which computes the difference as described in the equation above. Finally, we use plot_ramp to plot the median count rate from the left side, right side, and full frame image.

For our scattered light exposure, we see zodiacal light at a level of ~0.9e/s in later reads, with the scattered light component affecting the first several reads where the median count rate for the left side (orange triangles) is larger than the right side (green triangles). We can visualize this in 2 dimensions in the panel plot below, using plot_ima_difference_subplots.

In the panel plot, we see that sources (small galaxies) are visible in the difference images using this new method. Note that this may complicate the analysis of the spatial background (e.g. left versus right) for images with extended targets, such as large galaxies. In this case, users may wish to adjust the regions of the detector used for the ramp plots. We therefore recommend inspecting both the panel plots as well as the ramp fits for diagnosing any issues with the data.

fig = plt.figure(figsize = (20, 10))

ima_filepath = f'orig/{file_id}_ima.fits'


lhs_region = {"x0":50,"x1":250,"y0":100,"y1":900}
rhs_region = {"x0":700,"x1":900,"y0":100,"y1":900}

diff.plot_ima_difference_subplots(ima_filepath, difference_method='instantaneous', lhs_region=lhs_region, rhs_region=rhs_region)

<Figure size 2000x1000 with 0 Axes>

In this figure, we see that the ratio of instantaneous rate for the left versus right side of the image is ~1.0 for all but the first few reads (which are affected by scattered light). We choose to exclude reads with a ratio greater than 1.1 from the final science data, as they are the reads most strongly affected by scattered light. While this reduces the total exposure from 1403 to 1000 seconds, it removes the spatial component from the sky background.

4. Querying CRDS for the Reference File

In this notebook, we will run calwf3 to produce an IMA file, which we will use to remove bad reads. Before running calwf3, we need to set some environment variables.

We will point to a subdirectory called crds_cache using the IREF environment variable, which is used for WFC3 reference files. Other instruments use other variables, e.g., JREF for ACS.

os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_PATH'] = './crds_cache'
os.environ['iref'] = './crds_cache/references/hst/wfc3/'

The code block below will query CRDS for the best reference files currently available for our dataset and update the header keywords to point to these new files. We will use os to run terminal commands. In the terminal, the line would be:

crds bestrefs --files [filename] --sync-references=1 --update-bestrefs

...where 'filename' is the name of your fits file.

Warning: this cell may take a few minutes to complete.

raw_filepath = f'{file_id}_raw.fits'


print( f"Querying CRDS for the reference file associated with {raw_filepath}.")
command_line_input = 'crds bestrefs --files {:} --sync-references=1 --update-bestrefs'.format(raw_file)
os.system(command_line_input);

5. Correcting the Data

To address anomalies such as scattered light, satellite trails, and more in WFC3/IR images we can remove the individual affected reads and correct the science and error arrays accordingly. For our example image, we choose to exclude reads where the ratio of background signal is greater than 1.1 e-/s (see the notebook WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023) for a more complete demonstration of how we find this ratio).

5.1 Excising "Bad" Reads

In this section, we walk through the key steps to removing bad reads.

We select our excluded reads (in this case reads reads SCI,15 through SCI, 11) and use the "remove_reads" function to effectively subtract out the reads affected by TVB.

The key steps to this process are:

1. Run calwf3 on the RAW image file to produce an IMA and an FLT.
2. Take the IMA file and subtract the science array from the "bad" reads from the cumulative exposure.
   • Correct the final integration time accordingly.
3. Correct the error image extension to reflect the changes to the science array.
4. Update the final FLT headers, science array, and error array.

We walk through this process below.

STEP 1: Run calwf3

bad_reads = [11,12,13,14,15] #note, we leave the "zero read" (read number 16)

#### Remove any existing ima or flt products from working directory or calwf3 will not run 
for ext in ['flt','ima']:
    if os.path.exists(raw_filepath.replace('raw', ext)):
        os.remove(raw_filepath.replace('raw', ext))

 ##run CALWF3
calwf3(raw_filepath)

STEP 2: Correct the IMA product

Here we walk through the process of subtracting the reads affected by Earth-limb scattered light from the cumulative science extension of the IMA. While an FLT product is also created in step 1, we do not use this file due to the poor-quality ramp fit in the presence of TVB.

Overall, the steps are as follows:

1. Load the science data as a 3D array and the integration time as a 1D array from the science and time extensions of the IMA file. In addition, load the associated dark reference file.

ima_filepath_new = f'{file_id}_ima.fits'
ima = fits.open(ima_filepath_new)

cube, integ_time = diff.read_wfc3(ima_filepath_new)

dark_file = ima[0].header['DARKFILE'].replace('iref$', os.getenv('iref'))
dark_counts, dark_time = diff.read_wfc3(dark_file)

2. Convert science data from electrons/second to electrons.

cube_counts = np.zeros(np.shape(cube))
for img in range(cube.shape[2]):
    cube_counts[:,:,img] = cube[:,:,img]*integ_time[img]

3. Compute the difference between reads (the cumulative difference, in this case) to isolate counts and integration time from bad reads.

cube_diff = np.diff(cube_counts, axis = 2)
dark_diff = np.diff(dark_counts, axis = 2)
dt = np.diff(integ_time)

4. For each bad read, subtract the counts (and time) from ONLY that read from the total data.

final_dark = dark_counts[:,:,-1] 

final_sci = cube_counts[:,:,-1]
final_time = np.ones((1024, 1024))*integ_time[-1]


if (len(bad_reads) > 0):
    for read in bad_reads:
        index = cube_counts.shape[2]-read-1
        final_sci -= cube_diff[:,:,index]
        final_dark -= dark_diff[:,:,index]
        final_time -= dt[index]
        

STEP 3: Correct the Error Image

The errors associated with the raw data are estimated according to the noise model for the detector which currently includes a simple combination of detector readnoise and poisson noise from the pixel and which is intended to replicate the calwf3 error array calculation. Currently, the inital detector noise model (in electrons) is as follows, where RN is the readnoise.

$$ SCI_{1} = signal= (flux*flat + dark)$$$$ \sigma_{1} = \sqrt{(RN^{2}) + (flux*flat + dark)} $$

Note that our signal equation includes flat and dark corrections to revert the final science array (flux) produced by calwf3 to the original signal (in electrons) recorded in the detector.

The ERR array continues to be updated as the SCI array is processed. The dark error term is added in quadrature as the dark current is subtracted from the science array:

$$ SCI_{2} = signal - dark $$$$ \sigma_{2} = \sqrt{\sigma_{1}^{2}+\sigma_{dark}^{2}} $$

Next, the flat images are divided out of the science image. There are two flat-fields (p-flat and delta-flat), which are each divided out of the science image. The flat errors are combined using the correct error propagation method.

$$ flat = pflat * dflat $$

P-Flat (pflat) file:

$$ SCI_{3} = SCI_{2}/pflat $$$$ \frac{\sigma_{3}}{SCI_{3}} = \sqrt{ (\frac{\sigma_{2}}{SCI_{2}})^{2} +(\frac{\sigma_{pflat}}{pflat})^{2}} $$

Delta-Flat (dflat) file:

$$ SCI_{4} = SCI_{3}/dflat $$$$ \frac{\sigma_{4}}{SCI_{4}} = \sqrt{ (\frac{\sigma_{3}}{SCI_{3}})^{2} +(\frac{\sigma_{dflat}}{dflat})^{2}} = \sqrt{ (\frac{\sigma_{2}}{SCI_{2}})^{2} +(\frac{\sigma_{pflat}}{pflat})^{2} +(\frac{\sigma_{dflat}}{dflat})^{2}}$$

and finally it is converted to an error in electrons/second as:

$$ \sigma_{e^{-}/s}= \frac{\sigma_{4}}{exptime}$$

Evaluating this equation, we come up with an equation for error array calculation (in electrons/second) which we use in our Python calculations as follows:

$$ \sigma_{e^{-}/s}= \frac{\sqrt{ RN^{2} + (flux*(pflat*dflat) + dark) + \sigma_{dark}^{2} + (\sigma_{pflat}*flux*dflat)^{2} + (\sigma_{dflat}*flux*pflat)^{2} }}{pflat*dflat*exptime}$$

We do not address the error propagation through the nonlinearity correction step of calwf3 (NLINCORR), as it is complicated and sufficiently small to be outside of the scope of this notebook.

##getting variance terms
#### Readnoise in 4 amps
RN = np.zeros((1024,1024))
RN[512: ,0:512] += ima[0].header['READNSEA']
RN[0:512,0:512] += ima[0].header['READNSEB']
RN[0:512, 512:] += ima[0].header['READNSEC']
RN[512: , 512:] += ima[0].header['READNSED']


#### Gain in 4 amps
gain_2D = np.zeros((1024,1024))
gain_2D[512: ,0:512] += ima[0].header['ATODGNA']
gain_2D[0:512,0:512] += ima[0].header['ATODGNB']
gain_2D[0:512, 512:] += ima[0].header['ATODGNC']
gain_2D[512: , 512:] += ima[0].header['ATODGND']

#dark image
dark_im = fits.open(dark_file)

##And grabbing the flats
pflat_file = ima[0].header['PFLTFILE'].replace('iref$', os.getenv('iref'))
pflat_im = fits.open(pflat_file)
pflat = pflat_im[1].data

dflat_file = ima[0].header['DFLTFILE'].replace('iref$', os.getenv('iref'))
dflat_im = fits.open(dflat_file)
dflat = dflat_im[1].data


#Computing the final error

#### Variance terms

## poisson error: sqrt(signal), flux = final_sci, dark is converted to electrons
signal = (final_sci*(pflat*dflat)+final_dark*gain_2D)
## flat errors
dflat_err=dflat_im['ERR'].data
pflat_err=pflat_im['ERR'].data
#dark error
dark_err = (dark_im['ERR'].data)*gain_2D 
##final
final_err = np.sqrt(RN**2 + signal + (dark_err**2) + (pflat_err*final_sci*dflat)**2\
                  + (dflat_err*final_sci*pflat)**2)
    
final_err /= (pflat*dflat*final_time)
final_err[np.isnan(final_err)] = 0

#and finally the final science image (converted back to count rate)
final_sci /= final_time      

STEP 4: Update the FLT Product

Finally, we update the FLT product created in step 1 with the new science array, error array, and header keywords (such as the new total exposure time).

flt_filepath_new = f'{file_id}_flt.fits'
flt = fits.open(flt_filepath_new, mode='update')
 
#Updating the flt data

flt['SCI'].data = final_sci[5:-5,5:-5]
flt['ERR'].data = final_err[5:-5,5:-5]


# Updating the FLT header
flt[0].header['IMA2FLT'] = (1, 'FLT extracted from IMA file')
flt[0].header['EXPTIME'] = np.max(final_time)
flt[0].header['NPOP'] = (len(bad_reads), 'Number of reads popped from the sequence')

flt.flush()

Putting It All Together

Below, we have included the entire correction process (steps 1-4 above) in a function, for convenience.

def remove_reads(raw_filepath, bad_reads = []):
    '''
    From the final IMA read, subtract data of reads affected by anomalies.
    Compute a corrected science image and error image, and place them in the FLT product. 
    
    Parameters
    ----------
    raw_filepath: str
        Path to a RAW full-frame IR image fits file.
        
    bad_reads: list of int
        List of the IMSET numbers (science extension numbers) of reads affected by anomalies.
    '''
    
    flt_filepath = raw_filepath.replace('raw', 'flt')
    ima_filepath = raw_filepath.replace('raw', 'ima')
    
    #### Remove existing products or calwf3 will die
    for filepath in [flt_filepath, ima_filepath]:
        if os.path.exists(filepath):
            os.remove(filepath)
    
    #### Run calwf3
    calwf3(raw_filepath)
    
    #take calwf3 produced IMA, get rid of bad reads 
    ima = fits.open(ima_filepath)

    cube, integ_time = diff.read_wfc3(ima_filepath)

    dark_file = ima[0].header['DARKFILE'].replace('iref$', os.getenv('iref'))
    dark_counts, dark_time = diff.read_wfc3(dark_file)


    cube_counts = np.zeros(np.shape(cube))
    for img in range(cube.shape[2]):
        cube_counts[:,:,img] = cube[:,:,img]*integ_time[img]

    cube_diff = np.diff(cube_counts, axis = 2)
    dark_diff = np.diff(dark_counts, axis = 2)
    dt = np.diff(integ_time)

    final_dark = dark_counts[:,:,-1] 

    final_sci = cube_counts[:,:,-1]
    final_time = np.ones((1024, 1024))*integ_time[-1]


    if (len(bad_reads) > 0):
        for read in bad_reads:
            index = cube_counts.shape[2]-read-1
            final_sci -= cube_diff[:,:,index]
            final_dark -= dark_diff[:,:,index]
            final_time -= dt[index]



    ##getting variance terms
    #### Readnoise in 4 amps
    RN = np.zeros((1024,1024))
    RN[512: ,0:512] += ima[0].header['READNSEA']
    RN[0:512,0:512] += ima[0].header['READNSEB']
    RN[0:512, 512:] += ima[0].header['READNSEC']
    RN[512: , 512:] += ima[0].header['READNSED']


    #### Gain in 4 amps
    gain_2D = np.zeros((1024,1024))
    gain_2D[512: ,0:512] += ima[0].header['ATODGNA']
    gain_2D[0:512,0:512] += ima[0].header['ATODGNB']
    gain_2D[0:512, 512:] += ima[0].header['ATODGNC']
    gain_2D[512: , 512:] += ima[0].header['ATODGND']

    #dark image
    dark_im = fits.open(dark_file)

    ##And grabbing the flats
    pflat_file = ima[0].header['PFLTFILE'].replace('iref$', os.getenv('iref'))
    pflat_im = fits.open(pflat_file)
    pflat = pflat_im[1].data

    dflat_file = ima[0].header['DFLTFILE'].replace('iref$', os.getenv('iref'))
    dflat_im = fits.open(dflat_file)
    dflat = dflat_im[1].data


    #Computing the final error

    #### Variance terms

    ## poisson error: sqrt(signal), flux = final_sci, dark is converted to electrons
    signal = (final_sci*(pflat*dflat)+final_dark*gain_2D)
    ## flat errors
    dflat_err=dflat_im['ERR'].data
    pflat_err=pflat_im['ERR'].data
    #dark error
    dark_err = (dark_im['ERR'].data)*gain_2D 
    ##final
    final_err = np.sqrt(RN**2 + signal + (dark_err**2) + (pflat_err*final_sci*dflat)**2\
                      + (dflat_err*final_sci*pflat)**2)

    final_err /= (pflat*dflat*final_time)
    final_err[np.isnan(final_err)] = 0

    #and finally the final science image (converted back to count rate)
    final_sci /= final_time  


    #Updating the flt data
    
    flt = fits.open(flt_filepath, mode='update')

    flt['SCI'].data = final_sci[5:-5,5:-5]
    flt['ERR'].data = final_err[5:-5,5:-5]


    # Updating the FLT header
    flt[0].header['IMA2FLT'] = (1, 'FLT extracted from IMA file')
    flt[0].header['EXPTIME'] = np.max(final_time)
    flt[0].header['NPOP'] = (len(bad_reads), 'Number of reads popped from the sequence')

    flt.flush()

5.2 Comparing FLT Products

You may want to move your reprocessed images to a new directory, especially if you plan to run the notebook again.

reprocessed_flt = f'{file_id}_flt.fits'

original_flt = f'orig/{file_id}_flt.fits'

Now, we can compare our original and reprocessed FLT products.

image_new = fits.open(reprocessed_flt)
image_old = fits.open(original_flt)

fig = plt.figure(figsize = (20, 7))
fig
rows = 1
columns = 2

#add the total exptime in the title 
ax1 = fig.add_subplot(rows, columns, 1)
ax1.set_title("Reprocessed FLT image", fontsize = 20)
im1 = plt.imshow(image_new["SCI", 1].data, vmin = 0.8, vmax = 1.3, origin = 'lower', cmap = 'Greys_r')
ax1.tick_params(axis = 'both',labelsize = 10)
cbar1 = plt.colorbar(im1, ax = ax1)
cbar1.ax.tick_params(labelsize = 10)

ax2 = fig.add_subplot(rows, columns, 2)
ax2.set_title("Original FLT image", fontsize = 20)
im2=plt.imshow(image_old["SCI", 1].data, vmin = 0.8, vmax = 1.3, origin = 'lower',  cmap = 'Greys_r')
ax2.tick_params(axis = 'both', labelsize = 10)
cbar2 = plt.colorbar(im2, ax = ax2)
cbar2.ax.tick_params(labelsize = 10)

plt.rc('xtick', labelsize = 10) 
plt.rc('ytick', labelsize = 10) 

This new image was produced by subtracting the signal from the first 5 reads (not including the 0 read) from the final science data, reducing the effective exposure time from 1403 to 1000 seconds. While the total exposure is reduced from 1403 seconds to 1000 seconds (thus decreasing the overall S/N of the image), the background in the reprocessed image is now uniform over the entire field of view. We can see that the new FLT image is free of the Earth limb scattered light visible in the old FLT image.

print('The final exposure time after reprocessing is {}.'.format(image_new[0].header['EXPTIME']))

The final exposure time after reprocessing is 1000.00327.

6. Drizzling Nominal and Reprocessed FLT Products

In our example we use an exposure (icqtbbbxq) from image association ICQTBB020 acquired in visit BB of program 14037. This visit consists of two orbits of two exposures each, and we now download the three other FLTs in the visit (icqtbbc0q_flt.fits, icqtbbbrq_flt.fits, icqtbbbtq_flt.fits) and the pipeline drizzled DRZ product.

To produce a clean DRZ image (without cosmic rays), we can drizzle the four FLTs together (from the nominal exposures and the reprocessed exposure).

data_list = Observations.query_criteria(obs_id=OBS_ID)

Observations.download_products(data_list['obsid'], project='CALWF3',  
                mrp_only=False, productSubGroupDescription=['FLT','DRZ'])

Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/icqtbb020_drz.fits to ./mastDownload/HST/icqtbb020/icqtbb020_drz.fits ... [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/icqtbbbrq_flt.fits to ./mastDownload/HST/icqtbbbrq/icqtbbbrq_flt.fits ... [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/icqtbbbtq_flt.fits to ./mastDownload/HST/icqtbbbtq/icqtbbbtq_flt.fits ... [Done]
INFO: Found cached file ./mastDownload/HST/icqtbbbxq/icqtbbbxq_flt.fits with expected size 16583040. [astroquery.query]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/icqtbbc0q_flt.fits to ./mastDownload/HST/icqtbbc0q/icqtbbc0q_flt.fits ... [Done]

Table length=5

Local Path	Status	Message	URL
str47	str8	object	object
./mastDownload/HST/icqtbb020/icqtbb020_drz.fits	COMPLETE	None	None
./mastDownload/HST/icqtbbbrq/icqtbbbrq_flt.fits	COMPLETE	None	None
./mastDownload/HST/icqtbbbtq/icqtbbbtq_flt.fits	COMPLETE	None	None
./mastDownload/HST/icqtbbbxq/icqtbbbxq_flt.fits	COMPLETE	None	None
./mastDownload/HST/icqtbbc0q/icqtbbc0q_flt.fits	COMPLETE	None	None

nominal_file_ids = ["icqtbbc0q", "icqtbbbrq", "icqtbbbtq"]

nominal_list = []
for nominal_file_id in nominal_file_ids:
    shutil.copy(f'mastDownload/HST/{nominal_file_id}/{nominal_file_id}_flt.fits', f'{nominal_file_id}_flt.fits')
    nominal_list.append(f'{nominal_file_id}_flt.fits')
print(nominal_list)

['icqtbbc0q_flt.fits', 'icqtbbbrq_flt.fits', 'icqtbbbtq_flt.fits']

Next, we update the image World Coordinate System of the reprocessed image in preparation for drizzling.

updatewcs.updatewcs(reprocessed_flt, use_db=True)

Finally, we combine the four FLT images with AstroDrizzle.

astrodrizzle.AstroDrizzle('*flt.fits', output='f140w', mdriztab=True, preserve=False, build=False, context=False)

Comparing the new DRZ image made with the reprocessed FLT product against the original pipeline DRZ image, we see that the new DRZ image no longer includes scattered light but has a slightly lower S/N due the reduced total exposure time from 1403 to 1000 seconds.

DRZ_image = fits.getdata("f140w_drz.fits")
Orig_DRZ = fits.getdata('mastDownload/HST/icqtbb020/icqtbb020_drz.fits')

fig = plt.figure(figsize=(20,10))
rows = 1
columns = 2


#add the total exptime in the title 
ax1 = fig.add_subplot(rows, columns, 1)
ax1.set_title("Reprocessed DRZ Image", fontsize=20)
vmin,vmax = zscale(Orig_DRZ)
im1 = plt.imshow(DRZ_image, vmin=vmin, vmax=vmax, origin='lower', cmap='Greys_r')
_= plt.colorbar()

ax2 = fig.add_subplot(rows, columns, 2)
ax2.set_title("Original Pipeline DRZ Image", fontsize=20)
im2 = plt.imshow(Orig_DRZ, vmin=vmin, vmax=vmax, origin='lower', cmap='Greys_r')
_= plt.colorbar()

7. Conclusions

Congratulations, you have completed the notebook.

You should now be familiar with how to reprocess an observation affected by Earth limb scattered light by removing the affected reads from your science and error images.

Thank you for following along!

Additional Resources

Below are some additional resources that may be helpful. Please send any questions through the HST Help Desk.

• WFC3 Website
  • WFC3/UVIS CTE Website
• WFC3 Instrument Handbook
• WFC3 Data Handbook
• Instrument Science Report WFC3 2021-01

About this Notebook

Author: Anne O'Connor, Jennifer Mack, Annalisa Calamida, Harish Khandrika -- WFC3 Instrument

Updated On: 2023-05-16

Citations

If you use the following tools for published research, please cite the authors. Follow these links for more information about citing the tools:

• Citing numpy
• Citing astropy
• Citing astroquery
• Citing wfc3tools
• Citing drizzlepac

If you use this notebook, or information from the WFC3 Data Handbook, Instrument Handbook, or WFC3 ISRs for published research, please cite them:

• Citing the WFC3 Data Handbook
• Citing the WFC3 Instrument Handbook
• Instrument Science Report WFC3 2021-01
• Citing this notebook: Please cite the primary author and year, and hyperlink the notebook or WFC3 Library

---

Top of Page